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1. Introduction 

In many situations of interest, particles suspended in fluids cannot be modelled as simple 
point tracers. Both drops in gases and bubbles in liquids, and also solid powders in fluids, 
have a flnite size and their density is, generally speaking, different from the one of the 
advecting fluid. The description of their movement must then take into account the 
effects of inertia: this is why such objects are usually called inertial particles. Generally 
speaking, understanding the dynamics of these impurities [T[|2] is very relevant in several 
domains, ranging from geophysics [3l HI [5] to astrophysics [6l [7] , and from industry to 
biology [HIE]. 

Here, our attention will be focused on the sedimentation of inertial particles seeded 
in a given flow and subject to the action of gravity. Our main aim is to obtain an Eulerian 
description of the sedimentation process starting from the well-founded Lagrangian 
viewpoint for particle motion in the limit of small inertia [i.e. when collision events 
can be neglected). Although our main focus is the sedimentation process, our theory 
provides the whole detailed statistical information of particle motion. The probability 
density function of having a particle in a given position at a certain time is indeed 
available from our approach. Our theoretical machinery, which applies for a wide class 
of velocity flelds (either laminar or turbulent), is tested against numerical simulations 
(both direct numerical simulations and Lagrangian simulations) for a class of 2D cellular 
flows (static and time dependent). This speciflc example offers the possibility to analyse 
and discuss how the sedimentation process turns out to be extremely sensitive to the 
flow details. To be more speciflc, for the class of cellular flows, both an increase and 
a reduction of the rate of sedimentation with respect to the still-fluid case can occur, 
depending on the flow parameters like, e.g., flow geometry and pulsating frequency, as 
well as on the usual Stokes and Peclet numbers. Our flndings lead to the following 
conclusion of interest in the realm of applications {e.g. within environmental sciences 
in connection to the problem of aerosol dispersion and sedimentation): any attempt to 
model the effect of flows on particle sedimentation cannot avoid taking into account full 
details of the flow fleld. 

The paper is organized as follows. In section [2] the basic equations governing the 
time evolution of inertial particles under gravity and in a prescribed flow are given, 
together with the associated Fokker-Planck equation for the particle probability density 
function. In section [3] the problem of flnding the terminal velocity of the sedimenting 
particles is formulated in a Hermitian form and tackled via a second-quantization 
formalism. The resulting set of auxiliary equations to obtain the terminal velocity are 
presented in the limit of small Stokes numbers. In section H] the auxiliary equations are 
solved via direct numerical simulations in two dimensions and the obtained results are 
discussed. Conclusions and perspectives follow in section [5l The appendix is devoted 
to some mathematical and computational details. 
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2. General equations 

Let us consider the motion of a single, small, rigid, spherical particle of radius h immersed 
in an incompressible dimensional flow u{x, t) (with d >2). Even if some of our results 
are more general, the flow field will be assumed either steady or periodic in time (with 
period T), and periodic in space {i.e. cellular, with period L) [101 [H]. We shall focus our 
attention on the so-called Stokes regime, in which the surrounding flow is differentiable 
on scales of the order of b; we shall also neglect the feedback of the particle on the flow. 
The motion is thus influenced by gravity, buoyancy and drag [12], to which Brownian 
noise should be added in order to take into account the thermal fluctuations of the fluid. 
Moreover, we neglect the Basset, Faxen, Oseen and Saffman corrections [12] and other 
possible effects due to rotationality, high relative velocity or non-sphericity [T3] . 

To write the equation for the particle trajectory X{t), it is customary to introduce 
the covelocity V = X — [3u{X{t),t). Here, (3 = 3p{/{pf + 2pp) is an adimensional 
coefficient, where pp and pf are the densities of the particle and of the fluid, respectively. 
According to the ratio between the two densities, /3 ranges from (pf ^ Pp: heavy 
particles, like drops in gases) to 3 (pf ^ pp: light particles, like bubbles in liquids), and 
becomes 1 when the two densities are equal (and inertial effects absent). Therefore, the 
covelocity differs from the slip velocity {X — u{X{t),t)) by a term (1 — i3)u{X{t),t), 
which vanishes only for neutral particles. The concept of slip velocity is probably 
more familiar in connection with the study of relative acceleration. In our case, the 
use of covelocity is due to the fact that it strongly simplifies the equations of motion, 
properly keeping into account the added-mass effect without any need to introduce time 
derivatives of the external flow. However, u{x, t) being completely known, from the 
knowledge of V one can immediately find the real velocity X. 

The Stokes (response) time is defined as r = b'^/3i>(3, with u the kinematic viscosity. 
Denoting with g, k, and r){t) the gravity acceleration, the particle diffusivity and the 
standard white noise, respectively, we have: 



The study can be carried on in the corresponding phase space {x,v,t). Let us 
consider the probability that, at time t, the particle is at location x with a covelocity 
V, and denote it by p{x,v,t). This satisfies the Fokker-Planck equation associated to 
the stochastic differential equation ([1]): 



The linear operator C inside the curly brackets of ([2]) carries a factor k/t"^ in front of the 
Laplacian with respect to the velocity, which is the highest-derivative term. Therefore, 




(1) 




. 
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we may expect some singularity in the limits of vanishing diffusivity or large Stokes 
time. 

The terminal settling velocity of the inertial particle [lll|15] is defined as the average 
particle velocity 



w 



; J^^j^ j J^^ [v + Pu{x,t)]p{x,v,t) , (3) 

where the temporal integration is clearly omitted in the presence of steady flows. In 
general, such renormalized terminal velocity can differ from the bare terminal velocity, 
i. e. the particle settling velocity in a still fluid, whose value can easily be found to be 

w* = {l-(3)gr. (4) 

The difference turns out to be: 

Aw = w — w* = J^^'^ J y*^^ u{x, t)p{x, V, t) . (5) 

In what follows, we shall only take into account flows with zero mean (J dt J dx u{x, t) = 
0) and possessing odd parity with respect to reflections in the vertical direction. Such 
flows are relevant to analyse the vertical component Aw of the terminal velocity 
discrepancy because no mean contribution is present, and every (eventual) nonzero 
result, originated from the component of p antisymmetric in the vertical coordinate, is 
to be interpreted as due to preferential concentration in areas of rising or falling fluid. 

Our plan is thus to solve 1^, Cp = 0, for the phase-space density and to plug the 
result in expression ([5]) to find Aw. 



3. Analytical investigation at small Stokes number 

We shall first focus on situations in which the Stokes time is much smaller than the 
flow typical time scale, so that the particle adapts its own velocity to the one of 
the underlying fluid very rapidly. More precisely, denoting with L and U the typical 
length scale and velocity of the fluid respectively, in this section we shall only deal 
with small Stokes number, St = t/{L/U). As we are taking into account also gravity 
and diffusivity (besides inertia), a thorough adimensionalization of the problem also 
requires the introduction of the Froude and Peclet numbers, defined as Fr = [// 
and Pe = LU/k, respectively. As a first result, we see immediately that the vertical 
component of the bare terminal velocity (jlj), assumed as positive if pointing downwards 
and written in units of U, rewrites as 

= (1 -/3)StFr-2 . (6) 

Let us adimensionalize x with L, t with L/U and u with U, and let us denote 
the new variables with the same letters as before. Equation ([2]) can now be written in 
adimensional form after introducing an appropriate variable for the rescaled covelocity: 
the correct one turns out to be y = ^/tJ2kv. With this notation, the operator acting 
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upon the density p{x,y,t) (appropriately normalized in the new variables) becomes: 



C 



St 



-1 



d_ 



2 dy^dy^ 
d ' 



St 



+ st° 



-1/2 




_d_ 



d_ 
dt 



_d_ 

dx, 



+ st^/^ 




Pel-/3^ d 







2 Fr^ "^dy, 

= - St-Uo + St-^/^i + StM2 + St^/^s , (7) 

where G = g/g is the unit vector pointing downwards, and the expressions of the 
operators Ai {i = 0, . . . ,3) are easily identified from It is worth noticing that Ao 
is an Ornstein-Uhlenbeck operator acting on covelocity coordinates only, which leads to 
a Hermitian reformulation of the problem in terms of a second-quantization formalism 
for the covelocity variable. Technical details on the analytical computation are left to 
the appendix: here we only recall the essential steps. By introducing a power-series 
expansion in St for the particle probability density, we obtain a chain of advection- 
diffusion equations (see (Q) for the auxihary quantities 'ijjf^{x,t) {n G N) defined in 
the appendix, in which lower-order quantities act as source terms for higher-order ones. 
Such equations can be solved sequentially, at least numerically, once the external flow 
is given. As a result, expression ([5]) for the variation of terminal velocity (in units of U) 
rewrites 

L/U 



Aw = f2 St'+™ fdt^fdx u{x, t)^pl^,^^{x, t) , 

m=0 



(8) 



where the integral on covelocity has already been performed (see flA.7l) ). Note that 
only the vertically-antisymmetric components of V'4+2m(^' hereafter denoted with the 
superscript (in order to distinguish them from the vertically-symmetric counterparts 
^^^), contribute to the spatial integral in ([8]) for the vertical component Aw. Therefore, 
in order to study the leading contribution in Aw, quadratic in St (the equations needed 
to compute O(St^) are not reported here), it is sufficient to consider the system 



Mip, 



(S(a) 
4 



= -{l-P)Fi-^G 



d_ 

^ dx, 



(9) 



dx^ dXfj, 

MiJ^ = =^ = const. = , 

where we imposed uniform initial conditions (such that ipf, is a normalization constant) 
and denoted the adimensionalized advection-diffusion operator with 



-Pe- 



dx ^dx ^ 



(10) 



From ([H]), it is easy to understand that both "ifj^ and ip2 do not depend on Fr, and 
the former is clearly also independent of (3. This means that ■i/'f behaves like (1 — /3), 
therefore is proportional to both (1 — /5)^ and Fr~^. 
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Let us now focus on this leading contribution at order St^ in Aw, corresponding 
to the term m = in ([8]). From the analytical point of view, such considerations imply 
that we can introduce a non-universal (in the sense that it depends on Pe and possibly 
on other properties of the surrounding flow, e.g. its time frequency or space geometry, 
but not on St, Fr or /3) function 

/^W= [dt^ f dxG.uJx,t) '^\!^''''^\, , (11) 

such that we can write: 

Aw = [(1 - l3)St/Fr]^AW + 0{St^) . (12) 

Equation f[T^ summarizes the following two main analytical results for the variation of 
the actual terminal velocity in the limit of small inertia, i) It is influenced by gravity in 
the same way as the bare terminal velocity. Indeed, both w* in ([B]), and the coefficient 
in square brackets beside AW in ffTTj) . are proportional to Fr~^, i.e. to g. This result 
is not trivial, because e.g. it does not hold, in general, at higher orders in St. ii) It 
has the opposite parity with respect to the bare terminal velocity. Indeed, w* scales 
with (1 — (3) while the above-mentioned coefficient with its square. In other words, if a 
specific flow under speciflc conditions is found to increase the falling of heavy particles, 
then the same flow under the same conditions {e.g. same St, which may imply different 
L/U because r depends on P) must simultaneously decrease the rising of light particles, 
and viceversa. 

The dependence on molecular diffusivity is more difficult to obtain analytically. In 
particular, as already pointed out, the large- Pe limit is singular, and nothing precise 
can be said about intermediate Pe. On the contrary, for small Pe, it can be shown 
rigorously that the behaviour of AW is O(Pe^) or higher. As a particular case, we 
will show in the next section that, for the stationary square 2D Gollub flow, very well- 
deflned asymptotics are found numerically for both small and large Pe. In any case, one 
should keep in mind that ours is a series expansion in powers of St, whose coefficients 
may become very large when varying the other parameters. In order for this to work, 
St should be assumed small enough so as to make the following terms negligible with 
respect to the previous ones. 

4. Numerical investigation: the cellular flow 

From the analytical point of view, the results reported in the previous section are the 
only available, and all refer to the O(St^) correction. Proceeding at higher orders is not 
an easy task: e.g., the equation in (Q for ipl^^'^ would be much more complicated, and 
involves as a source term also the vertically-symmetric component tp'^^^^ , whose equation 
is in turn by far more cumbersome than the one for the corresponding antisymmetric 
part. Clearly, this problem also reflects on the numerical approach: solving such 
equations by Direct Numerical Simulation (DNS) is beyond our possibilities. The DNS 
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St Pe 

Figure 1. Left: comparison between LS (solid line) of ([T]) and the parabola (dashed 
line) resulting from the DNS solution of ([9]), for the variation of falling velocity as 
a function of the Stokes number at Pe = 5 and cj = 0. Right: dependence on the 
Peclet number for the variation of falhng velocity (solid hne) at O(St^); the dashed 
line represents the cubic asymptote found analytically, while the dotted line is the 
linear asymptote with a multiplicative coefficient derived from a fit. 



approach can however be apphed to ^/'g ^4 if the surrounding flow is simple 

enough. 

An accurate study has been performed for the 2D Gollub flow, or, more precisely, 
for a simple generalization thereof: 

ui = sin(/cxi) cos[x2 + sin(co't)] 
U2 = —k cos{kxi) sin[x2 + sin(co't)] , 

where the parameters k and u represent respectively the cell aspect ratio (vertical 
extension over horizontal one) and the pulsation of the oscillations in the vertical 
direction X2, assumed as positive downwards. The figures in this section all refer to 
the 2D Gollub flow in square cells {k = 1) for heavy particles {/S = 0) at the specific 
value Fr = 1. 

First of all, we test our analytical result on the O(St^) correction in the steady 
case (a; = 0). To do this, we performed a Lagrangian Simulation (LS) of the particle 
dynamics ([H) at Pe = 5 (see, e.g., [16], [TTj HH]). Additionally, we solved the system ([9]) 
for ifj^^^^ (and for ^/'f parallel) by means of DNS. In the left part of figured! for small 
(and intermediate) values of St, we report such LS and the parabola St^Al^, where 
AW is computed through (ITT]) with the ip^^^'^ found from DNS: the agreement is very 
good at least until St = 0.15. Moving at larger St, it is evident how the increase in the 
falling velocity reaches a maximum value at St ~ 0.4, then there is a crossover and a 
maximum reduction in settling takes place at St ^ 1. The asymptotic vanishing at large 
St is a general feature in the ballistic limit, in which the particle reacts so slowly to the 
variations of the surrounding flow that the latter has basically no influence on settling. 
Such limit is singular, as already pointed out, because a particle with infinite inertia 
in a finite-correlated flow is equivalent to a particle with finite inertia in a 5-correlated 
flow [19], [20] , thus a technique like uniform-coloured- noise approximation should be used 
to deal with this situation properly. On the contrary, we think that our perturbative 
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Figure 2. Left: same as in figure [T] (left) but with uj = 0.6; notice the decrease in the 
terminal velocity opposed to the static-carrier- flow case. Right: dependence on the 
pulsation uj for the variation of falling velocity (solid line) at O(St^), with Pe = 5. 



approach at small St should be able to capture at least the maximum in the plot, taking 
into account the following orders like if)^^^ or higher. Of course, this behaviour could 
not be caught correctly by our previous quadratic approximation, which is not able to 
describe inflection points and subsequent changes in convexity. 

Let us now turn to the dependence of AW on Pe. As already mentioned in the 
previous section, two well-defined asymptotics can be found, as shown in the right 
part of figure [H In particular, at small Pe, it can be shown by means of perturbative 
analysis that AW approaches f{k) Pe^ where f{k) = {k'^ + 3k^ - 2)/32{k'^ + if {i.e., 
f{k) > 4^ k > (vT7 — 3)/2 ~ 0.56). Note that this result is more stringent than the 
general behaviour O(Pe^) (or higher) reported in the previous section, and only applies 
to this particular instance of flow. On the contrary, the linear asymptote 0.19 Pe at large 
Pe comes from a fitting procedure, because this limit is singular. Notice, however, that 
this does not necessary imply that the terminal velocity diverges at infinite Pe: what we 
are discussing in this paragraph is only the behaviour of the O(St^) term, and clearly 
there is the possibility of a regularizing resummation if one also considers higher-order 
terms. 

In figure [2] we show that for a nonzero pulsating frequency the effect of the carrier 
flow may be the opposite (left panel). Indeed (right panel), the effective sedimenting 
velocity is decreased at small inertia in a window of values of uj. For larger frequencies 
{uj > 1.4) an increase in the sedimenting velocity occurs again. This clearly shows how 
sensitive to the flow details the terminal velocity may be. 

5. Conclusions and perspectives 

We have investigated how the terminal velocity of sedimenting particles is influenced 
by a background flow. In the limit of small St, we have reduced this problem to the 
solution of two coupled forced advect ion-diffusion equations. Such equations hold for 
generic velocities, laminar or turbulent, possessing odd parity with respect to reflections 
in the vertical direction. From our analysis the complete knowledge of the particle 
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motion is available: the probability density function of having a particle in a given 
position at a given time is directly accessible once the differential equations are solved. 
The determination of the terminal velocity is a particular product of our analysis. Apart 
from the limiting case of small Peclet numbers, one has to resort to numerical simulations 
to find the latter. This task is easily accomplished for two-dimensional laminar flows, 
which clearly lead to the following conclusions. The value of St alone is not sufficient to 
argue if the sedimentation is faster or slower with respect to what happens in still fluid. 
In particular, for square, static cellular flows, at both intermediate and high Pe, we 
found an increase of the falling velocity at relatively small St, and a reduction starting 
from St larger than some critical value. A similar transition is also found for a fixed St 
and assuming the cell to be oscillating in the vertical direction with frequency uj. We 
found the remarkable result that, appropriately tuning cj, a decrease in settling can also 
be obtained at small St, accompanied by an increase now occurring at sufficiently large 
St. We also studied the effect of the flow spatial geometry on the terminal velocity. 
In the limit of small Pe, the expression for the terminal velocity can be extracted 
analytically from our equations and its explicit dependence on the cell aspect ratio 
isolated. Changing the cell aspect ratio appropriately can cause a reduction in falling 
at small St. 

It would be interesting to make a comparison between our method and the well- 
known "continuum" approximation [Ij to obtain the particle probability density. In the 
latter case, diffusivity is neglected from the beginning. Some preliminary results indicate 
that at the lowest order in St the correction to the bare falling velocity coincide when 
computed with the two strategies. Substantial differences start to emerge at higher 
orders in St, where the "continuum" approximation probably ceases to hold, due to the 
formation of caustic patterns and to the phenomenon of clustering. 
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Appendix A. Mathematical and computational details 

In this section we provide some details of the calculation which leads to ([9]). 

A Hermitian reformulation of the problem is convenient. Let us define P]{y) = exp(— y^), 

the Gaussian kernel of ^o- -^oPt ~ ^- After decomposing the particle probability density 

as 

p(x, y, t) = p]^^{y)'^{x, y, t) , (A.l) 

1 /2 

our aim is to find the operators Bi such that AiP = — P| Biip (i = 0, . . . , 3). 

This can be achieved by means of a second- quantization algorithm. One can indeed 
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introduce the operators of creation and annihilation 

± 9 

in terms of which 

1 /2 

and the vacuum state turns out to be a~|0) = |0) — . 
Defining 

/Pe 1 d 



1 Q_ . _ /P^i-/3^ 

v^9x, ' V 2 Fr^ "^'^ 



the operators we are looking for read: 



_ d d 



Note that Bo can be interpreted as the occupation number, and the following 
commutation relations hold: 

[Bo, a J] = ±aj , [a;, a+] = 25^^ , [aj, a;^] = . 

On the other hand, B2 does not carry any creation or annihilation operator (A2 is 
independent of y) and vanishes when under investigation are very heavy particles in a 
stationary state. 

The following step consists in the development 

00 

i^{x,y,t) = Y,St''/'M^,y,t) (A.2) 

n=0 

and in the study of the set of equations arising at each (integer and half-i.) order in St: 

for n = 

I Biipri-i forn = 1 . 

Bilpn-l + B2llJn-2 for 71 = 2 

Bilpn-l + B2'll^n-2 + B^ljjn-S forn > 3 . 

Such relations can be solved recursively by exploiting the simple inversion formula for 
a generic function ^: 

1 

Therefore, one obtains 

Mx, y, t) = ^l{x, tm + Vn' {x, i)a+ |0) + . . . 

+ ^C^-""(^,^X---<J0), (A.4) 



i3o^ = <---<|o)^^' = T:<---<|o) 
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where 

rn'-^>'{x,t)= (A.5) 

+2k-^l,,^, {rn-r^'"''') for A; = 1, . . . , n - 2 
"MfcC-r""' for A; = n - 1, n 

and (. . .) implies a symmetrization on the repeated index {i.e. the sum of the possible 
permutations of the repeated index, divided by the number of such terms). Note that, by 
definition, the quantity ■ ■ -a^jO) equals the multidimensional Hermite polynomial 
of degree fc, Hj^{y). 

At each step one must also impose the corresponding solvability condition, which 
forbids the presence of states proportional to |0) on every right-hand side of expressions 
flA.Sp . in order to avoid inversion problems. These constraints give: 



27^<+i + i32C = VnGN. (A.6) 

By expressing each V'n+i; through ( 1A.5I) . as a function of a combination of ijjf^ for some 
m < n, expressions (1A.6P are to be interpreted as equations for the quantities ipl^. 



Together with the normalization condition J dx ijjf^ oc 6no, they can be solved recursively 
(numerically, or even analytically in some fortunate circumstances) once the basic flow 
u{x,t) is given. All such equations are of the advection-diffusion type and are forced 
by lower-order, already-solved quantities, with the exception of the equations for n = 
and n = 1, which are homogeneous and thus imply i/jq = const, and ijjf = under 
our assumptions. Moreover, one should notice that even-order equations are forced 
only by even-order quantities, and the same happens with odd n's: this means that 
ip2n+i = Vn G N; therefore, our expansion ( ]A.2[) reduces to a series in integer powers of 
St only. In we reported the first few equations corresponding to even n's, separating 
the components symmetric and antisymmetric in the vertical direction. Note that the 
spatial gradient and the flow field are vertically antisymmetric, while the vertical unit 
vector G (reminiscent of gravity) is symmetric in this sense. 

Starting from the quantities ipri^'^ with even n, we can reconstruct the variation of 



the adimensionalized terminal velocity by substituting back (1A.4P into (1A.2P into (lA.ip 



into ([5]), and by making use of the orthonormalization of the Hermite polynomials Hj^{y) 
in the integral in the y variable with weig ht pj^^(y) (note that 1 = H^): 

Aw = j dt^^ j dx j dyu{x,t)p{x,y,t) 

= j dt^^ j dxu{x,t) j dyp]''^{y)^{x,y,t) 
= jdt^^jdx uix, t) j dy e-y'l^ St"/'V^n(^, 2/> t) 



St"/2 I dt ^ I dx u{x, t) fdye 

n=0 J J J 



n=0 



X 
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k\ 

k=0 



E S*"^' E ^ /^^^ ^ ^^^n'-^" ^) X 

n=0 fc=0 '"^ 

X ldye-y"I^H^,{y) 
E St"^' E ^ ^ /d^ ^(^, 



n=0 fc=0 
oo 



m=2 

which corresponds to 



^ St"/2 fdt^fdx u{x, t)4jl{x, t) 
f^St"^ fdt^f dx u{x, t)^Ux, t) 

f2 St" fdt^fdx u{x, t)4;^ix, t) 

f^St- jdt^ ldxu{x,t)4'^\x,t) , (A.7) 
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